Mean field approximation of two coupled populations of excitable units 
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The analysis on stability and bifurcations in the macroscopic dynamics exhibited by the system 
of two coupled large populations comprised of N stochastic excitable units each is performed by 
studying an approximate system, obtained by replacing each population with the corresponding 
mean-field model. In the exact system, one has the units within an ensemble communicating via the 
time-delayed linear couplings, whereas the inter-ensemble terms involve the nonlinear time-delayed 
interaction mediated by the appropriate global variables. The aim is to demonstrate that the bifur- 
cations affecting the stability of the stationary state of the original system, governed by a set of 4iV 
stochastic delay-differential equations for the microscopic dynamics, can accurately be reproduced 
by a flow containing just four deterministic delay-differential equations which describe the evolution 
of the mean-field based variables. In particular, the considered issues include determining the pa- 
rameter domains where the stationary state is stable, the scenarios for the onset and the time-delay 
induced suppression of the collective mode, as well as the parameter domains admitting bistability 
between the equilibrium and the oscillatory state. We show how analytically tractable bifurcations 
occurring in the approximate model can be used to identify the characteristic mechanisms by which 
the stationary state is destabilized under different system configurations, like those with symmetrical 
or asymmetrical inter-population couplings. 
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The onset and mutual adjustment of collective rhythms 
are regarded as the dynamical paradigm for the macro- 
scopic phenomena in a wide range of biological and in- 
organic systems. Such a framework has already proven 
indispensable for understanding the normal and patho- 
logical patterns of brain activity 0, [|J , coordination of 
cellular clocks governing the circadian rhythms Q, the 
mechanisms regulating heartbeat [|| or lying behind cer- 
tain forms of social behavior Q , entrainment of electro- 
chemical oscillators rah as well as the dynamics of Joseph- 
son junction circuits [7j and the arrays of coupled lasers 
Q . The emergence of macroscopic rhythms in ensembles 
of oscillating units is mediated by the synchronization 
based self-organization !§] . The latter is often influenced 
or facilitated by noise on one hand , while on the 

other hand, the interaction over the appropriate commu- 
nication channels is typically susceptible to transmission 
delays or there may be a time lag due to the system 
components' latency in response to input variations [l2j]. 
A pervasive idea in nonlinear dynamics is to treat an 
assembly exhibiting a collective mode as a macroscopic 
oscillator [la ], which could in turn be subjected to an 
external drive or be exposed to a single or multiple col- 
lective rhythms from other populations. In this context, 
an important issue is to consider the relationship between 
the ensembles's global variable and the external forcing 
or that between the corresponding global variables. 
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In terms of the dynamical complexity of the observed 
behavior and the methods available for the analytical 
study, one has to make a distinction between the cases 
where the populations are built of self-sustained (au- 
tonomous) oscillators or the excitable units. In the for- 
mer instance, it is possible to obtain a more compact 
description of the interacting ensembles' dynamics by ap- 
plying the phase reduction techniques |14T - [la |. Given 
that the phase cannot be attributed to the system re- 
siding at the equilibrium, excitable populations are not 
amenable to such methods. Nonetheless, on the level 
of elementary behavior associated with the macroscopic 
variables, populations containing the excitable or self- 
oscillating units undergo qualitatively similar forms of 
dynamics. In particular, the ensuing collective modes 
may synchronize [l7], HI], become phase- locked or get 
suppressed by the action of the coupling delay (delay- 
induced amplitude death) 19]. Beyond such simple 
cases, there are more complex forms of collective behavior 
tied exclusively to populations of interacting oscillators. 
A few prominent examples include the self-organized 
quasiperiodicity FH1 an d the partially synchronous chi- 
maera states [20L |21| , which have been found to emerge 
in systems of identical phase oscillators under the action 
of external forcing or by coupling to another population, 
respectively. The former regime is characterized by the 
frequency of the collective mode being distinct from that 
of the single elements, while the other involves a bro- 
ken symmetry between the dynamics of two interacting 
populations. 

In this study, the focus lies with the two delay-coupled 
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populations of identical excitable units modeled by the 
Fitzhugh-Nagumo elements. The behavior of the latter 
is representative for the type II excitability [22J , which in 
contrast to type I, lacks a sharp threshold in a sense that 
the amplitude of the response depends continuously on 
the size of the applied stimulus. Though the considered 
framework is quite general, the basic motivation admit- 
tedly draws from the observations on neuronal assem- 
blies, with the adopted model of local activity typically 
invoked in such a context. The analysis of the underlying 
system dynamics may be approached from two different 
angles. For one, a numerical study can be carried out 
to look for the states of the increasing dynamical com- 
plexity. Instead, we take on a strategy that consists in 
examining how well is the behavior of the exact system 
matched by that of the coupled mean-field (MF) sys- 
tems, having derived the MF model as an approxima- 
tion for the activity of a single ensemble. The concept 
aims to fully exploit the analogy between the assemblies 
and the macroscopic oscillators, such that the original set 
of equations for the microscopic dynamics is reduced to 
a flow which describes the evolution of the global vari- 
ables, incorporating the cross-population interaction in 
a natural way. An important ingredient for the setup 
is that both the intra- and the cross-population coupling 
terms include the transmission delays. Note that the lay- 
out with two populations may constitute a paradigm, or 
rather serve as a nucleus for the "network of networks" 
[l7l I2H [23[ | , which can be realized as a hierarchy of mul- 
tiple networks, or it could be thought of as an idealiza- 
tion for a single network with a strong modular struc- 
ture and a large number of elements in each community 
(subnetwork). Both configurations are common in bi- 
ological systems [l7|, ranging from the cellular level to 
the distributed anatomical areas of the brain, and also 
encompassing the populations of cells responsible for the 
rhythmic activity in heart, kidney, pancreas, to name but 
a few. As for the comparison with the MF model, the 
attempts at providing a reduced description instead of 
using the complete set of equations for each and every 
population constituent, have a particularly long history 
within the neuroscience (24T - I261 ]. Apart for the gains on 
the modeling side, they have initially been instigated by 
the finding that the EEG and MEG recordings may be 
linked to an average behavior, viz. the massively summed 
action potentials emitted within the strongly coupled, 
but remote cortical areas [HI, H3] ■ Though the given ap- 
proach inevitably includes simplifying assumptions that 
eventually constrain the repertoire of possible system be- 
haviors just to periodic motion, some of the realism may 
readily be sacrificed for a more parsimonious represen- 
tation if the emergence of the collective mode and the 
related dynamics are reproduced with sufficient fidelity. 

The mean field approximation has been applied on sys- 
tems of excitable units with noise but with no time-delay 
for example in Otherwise a type of MF 

approximarion was devised in |32| and [33| and applied on 
large clusters of noisy neurons with time-delayed interac- 



tion in 34]. However, the approximations made in these 
papers resulted in a system of equations that is still to 
large to be analyzed analytically, so that the approximate 
system must be studied numerically. We shall utilize an 
approximate system of only two DDDE, introduced in 
[3a |. for the dynamics of the mean fields for each of the 
two populations. Such a simple system allows analytical 
treatment of bifurcations and the parameter domains of 
stability of the stationary states which turn out to be in 
a quite good agrement with the exact complex system. 

The key set of issues addressed in this study amounts 
to identifying the conditions for the stability of the sta- 
tionary state, the onset of the collective mode, bistabil- 
ity between the equilibrium and the oscillation state, as 
well as the time-delay induced suppression of the collec- 
tive mode. One notes that the applied term "collective 
mode" here implies the existence of a limit cycle for the 
total system of interacting populations. Though the in- 
tention is not, or rather cannot be to account for any 
experimental observation of such phenomena, some ele- 
mentary comparison can still be drawn. For instance, the 
notion that the emergence and the synchronization prop- 
erties of collective rhythms arising in the macroscopic 
neural populations are critically influenced by the cou- 
pling strength and the interaction delay [27j has its clear 
analogue in the results we arrive at. Consistent with the 
stated objectives, the study of the approximate system is 
concerned with the local bifurcation analysis, carried out 
analytically and corroborated by the numerical means, 
to determine i) the parameter domains of stability of the 
steady states, ii) the scenarios for the emergence or the 
suppression of the collective mode, and Hi) the parame- 
ter domains admitting the bistability between the equi- 
librium and the oscillatory state. 

The paper is organized as follows. In Section HI the 
details of the exact model of interacting populations are 
laid out in parallel with the derivation of its MF coun- 
terpart. Section [XT] is focused on the local bifurcation 
analysis of the approximate model, providing for the an- 
alytical results. In Section IIII1 we demonstrate that the 
approximation based on two coupled MF systems is able 
to accurately predict the behavior of the exact system 
in terms of the stability of the equilibrium, as well as 
the onset and the suppression of the collective mode. It 
is also pointed out how different system configurations 
affect the scenarios for the emergence of the oscillatory 
state and influence the parameter domains supporting its 
coexistence with the equilibrium. The results are briefly 
summarized and discussed in the concluding section. 



I. BACKGROUND ON THE EXACT MODEL 
AND DERIVATION OF ITS MF COUNTERPART 

A. Details of the exact model 

Each population comprises a collection of N identical 
Fitzhugh-Nagumo elements [22I l37j , whose dynamics is 
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given by 

edx i: i = (xi,i - - y i: \ + h)dt + ^^-x 

N 

y")fcj,l(* - T in,l) - %i,l(t)]dt + 5c,l X 

arctan[X2(t — r c ,i) + 62]^, 
= (a;,-,! + 61)* + y/2EhdW itl 
edx l:2 = (x it 2 - x 3 l 2 /3 - y ia + Ja)d* + 

AT 

/Jfoj.aft - T in , 2) - x lt2 (t)]dt + g Ct2 x 
3=1 

arctan[_X"i (t — t c2 ) + bi]dt, 
dy it 2 = (Xi,2 + h)dt + y/2D^dW h2 , (1) 

where the subscripts k = 1,2 specify the population, in- 
dices i = 1, ..N denote a particular unit within the pop- 

N 

ulation, and Xk = (1/iV) ^2 x% t k stand for the macro- 

»=i 

scopic variables that typify the global population be- 
havior. The small parameter e = 0.01 imposes a wide 
separation between the characteristic time scales for the 
evolution of x^k and y^k- In the context of neuronal 
activity the set of fast variables embodies the mem- 
brane potentials, whereas the slow- variable set is sup- 
posed to account for the gross kinetics of the potassium 
ion-gating channels. In the absence of an external stim- 
ulation Ii = I2 = applies. The impact of a noisy 
background activity is reflected by the %/ 2DdWi terms, 
which represent the stochastic increments of the indepen- 
dent Wiener processes specified by the noise amplitude 
D, expectation values (dWi) = and the correlations 
that satisfy (dWidWj) = Sijdt for each population. 

Owing to the system configuration, the local dynam- 
ics involves two types of interactions, each characterized 
by the coupling strength and the delay. The respec- 
tive parameters associated with the intra-ensemble terms 
are gi n ,k and Ti n ,k, while the cross-population terms are 
awarded g c> k and t c ^. Within the populations, the ele- 
ments communicate via the simple linear (diffusive) cou- 
plings, such that Ti n may account for the transmission 
delays due to finite rate of signal propagation or the la- 
tency in unit responses. Given the objectives stated in 
the Introduction, it is not unjustified to make use of some 
simplifying assumptions, like the all-to-all pattern of in- 
terconnections and the uniformity of coupling strengths 
inside the ensembles, which are the abstractions often 
invoked in the relevant literature [38j]. As for the cross- 
population interactions, at the current stage no particu- 
lar model is considered to be preferred over the others. 
However, we make use of an analogy to neural systems 
by noting how a variety of models display a common fea- 
ture. Stated in the language of neuroscience, the evoked 
postsynaptic potentials can be expressed in a symbolical 
form h = s (£> m, where m refers to an average density 



of presynaptic input arriving from the transmitter pop- 
ulation, and s presents the threshold-like response of the 
neurons of the receiving population (2?J ■ Adhering to this 
concept, the output of the transmitter population is inte- 

N 

grated by the macroscopic variables X^ = (1/iV) x i,k, 

i=l 

which reflect the global behavior in a sense that the bet- 
ter the synchronization among the constituent elements, 
the larger the amplitudes of X^. In terms of the nonlin- 
ear threshold function, there is a degree of arbitrariness, 
so the arctan form applied here is as good a choice as any. 
Unlike the interactions within the populations, which are 
characterized by the specific strengths per link, the inter- 
population terms involve the cumulative strengths, con- 
sistent with the idea of viewing each population as a sin- 
gle macroscopic oscillator. The meaning of the parameter 
b is explained in more detail further below. The bidi- 
rectional couplings between the ensembles, being either 
symmetrical or asymmetrical, may be important from the 
aspect of neuroscience, given that the brain connectivity 
patterns are known to exhibit a large portion of recipro- 
cal interactions [261 ]. On the level of local dynamics, the 
parameter b plays a key role as it modulates the unit's 
excitability. For an isolated unit in the noiseless case, the 
condition \b\ = 1 determines the Hopf bifurcation thresh- 
old, above which the system possesses a unique equilib- 
rium, whereas below it one finds a limit cycle. Selecting b 
slightly above 1, like the value b — 1.05 held throughout 
the paper, the population elements are poised quite close 
to the Hopf threshold, which gives rise to an excitable 
behavior. In such a regime, an adequate stimulation, be 
it by the noise or the interaction terms, may evoke large 
transients within the fast variable subspace before the 
ensuing trajectories converge back to rest. Note that in 
the scenario where noise acts in the slow subsystem, the 
elicited limit cycles are just the precursors of the deter- 
ministic ones [39(. Turning back to the role of b in the 
cross-population coupling, it is seen to ensure that the 
largest contribution to the interaction term comes from 
the global states lying farthest away from the equilib- 
rium. 



B. Background and the formulation of the MF 
approximation 

Deriving the MF approximation, we aspire for a highly 
reduced set of nonlinear DDE instead of the original sys- 
tem (JXJ) comprised of a large set of nonlinear SDDE. 
Though a simplified representation, the MF model should 
still be able to reproduce with sufficient accuracy the lat- 
ter 's behavior regarding the stability of the steady states, 
the scenarios for the onset of the collective mode and its 
suppression under the action of the cross-population cou- 
pling delay. The MF treatment draws on the all-to-all 
type of connectivity among neurons within each popula- 
tion, incorpora ting the thermodynamic limit N — > 00 in 
a natural way [25|]. In order to build a MF model, two 
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different approaches are available to proceed with: one 
may either consider the time-dependence of a hierarchy 
of probability densities according to the Fokker-Planck 
formalism, or may focus on the evolution of cumulants, 
whereby the full density of states is factorized into a series 
of marginal densities. The latter alternative is preferred, 
as it allows for a number of convenient approximations 
to be introduced in a controlled fashion [25(- Note how 
one is bound to make some approximations for the non- 
linearity of the original system, given that the cumulants 
of the particular order are usually linked to those of the 
higher order, which apparently renders the underlying 
series unclosed. The way to resolve this issue consists 
in truncating the series by a form of a closure hypoth- 
esis. Such hypothesis typically integrates the cumulant 
approach with the Gaussian approximation pol l4lj , re- 
calling that the Gaussian distribution has vanishing cu- 
mulants above the second order. 



Confined to a single population, the Gaussian approx- 
imation involves two elementary prepositions: first, that 
the instantaneous distributions of local variables P(xi) 
and P(yi) are Gaussian, and second, that the ensemble 
averages at any given moment coincide with the expec- 
tation values of the appropriate distributions in a sense 

N N 

(i/N) E[p( Xi )}, (i/N) Ey«« E[p{ y% )} nam. 

»=1 i=l 
If the two stated conditions are met, all the cumulants 

above the second order are supposed to vanish. Let 
us briefly comment on the constraints which these con- 
ditions impose on the system parameters. As for the 
first point, the Gaussian distribution of local variables 
is maintained if the noise amplitude obeys D << 1. 
Nonetheless, the strong law of large numbers [42[ im- 
plies that the second condition concerning the ensemble 
averages is fulfilled exactly in the thermodynamic limit 
N —> oo if the involved stochastic processes are indepen- 
dent (gi n << 1). However, the numerical results pre- 
sented further on indicate that the MF approximation 
remains valid if the two latter conditions are relaxed, viz. 
when there is non-negligible interaction in the finite-size 
systems, provided that the requirement for not too large 
a noise amplitude is satisfied. 



In the following, we outline the key steps in the deriva- 
tion of the MF model for the activity of an interacting 
assembly. The derivation presents a slight generalization 
of the one presented in [35[ . To begin with, note that the 
cross-population coupling terms involve only the aver- 
age dynamics of the respective transmitter populations. 
This means that the focus should really lie with the inter- 
nal ensembles' dynamics, treating them temporarily as if 
they were independent, while subsequently including the 
inter-population interaction. Therefore, we confine fur- 
ther presentation to a single population, whose dynamics 



is extracted from ([T]) by setting g c ^\ or g c ^ to zero 

N 

edxi = (xi - xf/3 - y t )dt + -j^y~][xj(t - r in ) - Xi (t)]dt, 

3=1 

dy % = (xi + b)dt + V2DdWi, (2) 

Given that the distributions of the stochastic local vari- 
ables are assumed to take on the Gaussian form, one can 
fully characterize them by the set of the first and sec- 
ond order moments, which includes the mean values, the 
variances and the covariance. The mean values applied 
here 

N 

m x (t) = ( Xi (t)) = lim (l/N)Y y X l (t) 

»=1 
N 

m y (t) = Mt)) = lim (I/N) Vi(t) (3) 

i=l 

should strictly speaking be distinguished from the global 
variables X and Y considered earlier for the large, but 
still finite-size populations. The angled brackets are gen- 
erally used to denote averaging over the units making up 
the ensemble, whereas m x and m y are reserved solely for 
the averages of the local variables. Before introducing 
the second order moments, it is convenient to define the 
deviations from the mean n Xi (t) = (xi(t)} — Xi(t) and 
n Vi (t) — (yi(t)) — yi(t), which obey the Gaussian distri- 
butions and are independent between the single elements. 
Then the appropriate variances read 

s x (t) = {nl i (t))={({x i (t))-x i (t)) 2 ) 
«.(*) = «(*))=<«»(*)> -W(*)) 2 >. ( 4 ) 
whereas the covariance is given by 

u(t) = {n Xi {t)n Vi {t)) = (((xiity-XiitMyiW-yiQ))). 

(5) 

The evolution of the distributions' means m x and m y 
is obtained by performing the ensemble averages over 
the system ([2]), while the expressions for the dynamics 
of s x ,s y and u follow from explicitly taking the time 
derivatives of the definitions (j4j and ([5]). Note that 
the latter calculation also involves the derivatives of the 
compound functions of the stochastic variables such as 
d(x 2 )/dt and d(y 2 )/dt, where one is required to ap- 
ply the Ito's chain rule. As for the higher order aver- 
ages, like (xf) and (xf), it is necessary to tie them to 
the first and second order moments. In the simplest 
cases, this is accomplished by using the definitions (|4| 
and ([5]), while in most instances one arrives at the re- 
quired relations by setting the higher order cumulants 
!43j to zero, e.g. (xf) c = (zf) - 3(x 2 i )(x i ) + 2( Xl ) 3 = 0, 
{xhilc = (xfyi) ~ ( x i)(yi) - 2(x l )(x l y l )+2(x l ) 2 (y l ) = 0, 
and similar for (x^yi) c — and (xf) c — 0. The ensu- 
ing auxiliary formulas for the higher order averages then 
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read 



(&<) 


2 

= s x + m x 






= ml + 3m x s x 




(4) 


= m x + 6m x s x + 


34 




— u + m x m y 






= m y s x + m y m x 


+ 2m x u 




= 3s x u + 3m x u -\ 


- m y m x + 2>m x m y s x 



(6) 

After a series of steps which are too lengthy to convey 
in full detail, the closed system of equations for the first 
and second order moments finally becomes 



dm x (t) 
: dt 



= m x (t) - m x (t) 3 /3 - s x (t)m x (t) - m y (t)- 
hn{m x (t -T in ) - m x (t)) 

= m x (t) + b 

= s x (t)(l - m 2 (t) - s x (t) - g m ) - u(t) 
= u(t) + D 



dm y (t) 

dt 
£ ds x (t) 

2 dt 

1 ds y (t) 

2 dt 

S =^(1 - mt(t) - s x (t) ~ .9m) - - e s y (t) + s x (t). 

(7) 

Note that (|7|) comprises a set of deterministic delay equa- 
tions, where the impact of noise is absorbed into its am- 
plitude D. Recalling the Introduction, one of the objec- 
tives has been to carry out the bifurcation analysis on the 
MF model analytically. However, the system ([7]) is still 
sufficiently complex to defy such a treatment. To ensure 
that the bifurcation analysis is analytically tractable, we 
consider an additional approximation which concerns the 
relatively fast relaxation of the second order moments. 
Given that the characteristic time scales, at least for s x 
and u, are dominated by the small parameter e, one may 
substitute their full dynamics by the stationary values 
reached when s' x = 0, s' y = and u = are satisfied. 
Though a crude approximation, it is not an uncommon 
one [25|, |4l|. In the language of neuroscience, the net 
result it yields is comparable to translating the initial 
MF model into an effective neural- mass model [27|, the 
former (latter) associated with the system of five (two) 
equations. Nevertheless, whether this is justified or not 
strongly depends on the main objectives of the study, 
which here concern the stability of the stationary state, 
the onset of the collective mode and its suppression in 
an amplitude death-like phenomenon [TT| . As it stands, 
the described modification to the MF model should not 
substantially affect the latter set of issues, since the infor- 
mation supplied by the second order variables, like that 
on small fluctuations around the collective synchronous 
state, appears redundant in such a context. This is cor- 
roborated later on by the results indicating an agreement 



between the behaviors of the exact and the MF approxi- 
mation. 

To complete the MF approximation for the dynam- 
ics of the two interacting populations, one should take 
into account the inter-ensemble interactions initially left 
aside, arriving at the following set of four equations 

dm xA (t) m x i(t) 3 m xA (t) 
e = m Xt i(t) — (1 - g in ,i- 



dt 



m xA (t) 2 + yj (g tnA - 1 + m xA {t) 2 ) 2 + 4£>i)- 

TO y,lW + 9m,lK,l(t - Tin,l) - m x> i(t)) + 

g Ct i arctan(m X:2 (i - t c ,i) + b 2 ) 
m Xi2 (t) 3 m xA (t) 



dm v i(t) , , 

dm xfi {t) 



m.2 



(*) 



(1 - Sin, 2 



dt —.-w 3 2 

m x , 2 (i) 2 + sj (g m .2 - 1 + m x , 2 (t) 2 ) 2 + 4D 2 )- 

™ y ,2(t) + gin.i{rn x . 2 (t - r irl)2 ) - m x #(t))+ 
g c ,2 arctan(m :r4 (t - r ca ) + 61) 

±^ J - = m x>2 (t) + b 2 

dt 

(8) 

Note that for Di = D 2 = 0, the obtained system strongly 
resembles the case of two interacting Fitzhugh-Nagumo 
elements subjected to the delayed feedback. Another 
point is that the isolated populations (g c ,i — 3c, 2 = 0) 
can be shown to exhibit the excitable-like dynamics un- 
der the variation of D and r. By this is meant that apart 
from the small amplitude oscillations about the equilib- 
rium, there may also be large excursions of the global 
potential, this reflecting the crucial feature of the ex- 
act system. In our previous paper, it has already been 
demonstrated that the MF model of a single assembly 
is able to accurately predict the qualitative behavior of 
the exact system [35| . This refers to a sequence of local 
bifurcations under variation of D, Ti n and gi n , which can 
be used to highlight the parameter domains giving rise 
to oscillatory states or those that lead to the amplitude 
death 36]. In addition, the MF model of an isolated 
ensemble has also been found to reflect the global bifur- 
cation imminent to the onset of clustering in the exact 
system |44j |. 

Before proceeding to the main results, several brief re- 
marks on the applied numerical integration schemes are 
in order. The time series for both the exact and the 
approximate models are obtained by implementing the 
Euler method with the fixed time step At — 0.005 in the 
former, and At — 0.01 in the latter case, having verified 
that no changes occur for the smaller At. Also, on either 
occasion, we have adopted the standard and physically 
plausible initial functions, based on the assumption of 
the units evolving independently within the time interval 
t G [-T m i„,0], where T min = mm{r in ,i, T in<2 , r C) i, r C)2 }. 
This effectively amounts to integrating the systems ([1]) 
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and dS]) by disregarding any interaction for t € [— T min , 0], 
with the initial conditions in each instance taken in the 
vicinity of the fixed point. The results for the exact model 
refer to populations made up of N — 200 elements, but 
have been verified to persist if the larger assemblies are 
considered. 



II. ANALYTICAL RESULTS OF THE LOCAL 
BIFURCATION ANALYSIS OF THE 
APPROXIMATE SYSTEM 



In the two following sections, we first provide the de- 
tails of the local bifurcation analysis performed on the 
approximate system and then examine whether and how 
well do these results match the behavior of the exact sys- 
tem, whereby the latter dynamics is represented by the 
typical sample paths obtained from numerical integration 
of (HJ) for the sufficiently large N with D X ,D 2 ^ 0. On 
the first part, the analysis covers the stability of the at- 
tractor states for the total system of coupled populations, 
such that both of them are either found lying in the equi- 
librium or exhibiting oscillations. The main focus is on 
the stability of the fixed point and its destabilization un- 
der variation of the cross-population coupling strengths 
and delays. Apart for the onset of the oscillatory state, it 
is also considered how the coherent rhythms may become 
suppressed, this primarily attributed to the action of the 
inter-ensemble time lags. As a final matter, we demon- 
strate the existence of the parameter domains admitting 
the bistable regime, where the stationary and the oscil- 
latory state coexist. Altogether, an inference confirmed 
later on is that the MF approximation can capture the 
behavior of the exact system much better if the collective 
dynamics is such that the deterministic component, con- 
trolled by the coupling strength and time delay, prevails 
over the stochastic component. The points enumerated 
above exhaust the corpus of problems that may approx- 
imately be treated by the local bifurcation theory, in a 
sense of explaining the qualitative changes arising in the 
system's asymptotic dynamics due to parameter varia- 
tion. Outside the scope remain the more complex phe- 
nomena occurring for larger D-s and r-s, which could 
cause the behavior of single units within the populations 
to become substantially stratified. Such issues would fall 
under the notion of stochastic bifurcations [42| , meaning 
that one should consider how the parameter modifica- 
tion influences the changes of the respective stationary 
distributions of the local variables. 

Since we discuss the scenarios with symmetrical and 
asymmetrical cross-population couplings, as well as the 
setups where the inherent ensemble dynamics is the same 
or distinct, the analytical results of the local bifurcation 
analysis on the system of interacting MF models are pre- 
sented in most general terms with respect to the system 
parameters. First, it is established that the system © 



possesses a unique equilibrium given by 

b k 
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with k = 1,2. The local stability of (|9]) depends on 
the roots of the characteristic equation of the system 
([8]). To obtain the latter, one linearizes (JU) around 
the equilibrium, assuming that the deviations are of the 
form Sm Xik (t) = A k e xt ,Sm ytk (t) = B k e xt and 6m Xtk (t — 
T in,k) — A k e x ( t ~ Tin - k \ This results in a set of algebraic 
equations for the coefficients A k and B k , which has a 
nontrivial solution only if 



Ai(A)A 2 (A) - \ 2 g c , l9c .2e- x(T - 1+T ^ ] = 



(10) 



is fulfilled, where A fe (A) = -\F k + e\ 2 - g in , k \e- AT ^ k + 1 
with F k = F k (g inkl b k , D k ). The condition (fTU)) poses the 
desired characteristic equation, whose being transcenden- 
tal reflects the presence of (multiple) time delays in ([5]). 
Though (|10p has an infinite number of roots, it is well 
known how there may be only a finite number of excep- 
tional roots equal to zero or with a zero real part (45l . l4q ] . 
One recalls that tangent to the subspace spanned by the 
associated eigenvectors lies the center manifold 47], 13 , 
where the qualitative features of the system's dynamics, 
such as the local stability, are contingent on the nonlinear 
terms. 

Bifurcations of the stationary state take place for the 
parameter values where the roots of (fT0|) cross the imag- 
inary axes. Given that Eq. (fTU]) does not admit the 
possibility A — 0, we look for the pure imaginary roots 
of the form A = iuj, adopting uj to be real and positive. 
Substituting for A in (jTDJ) , one obtains 



[-ioj(Fi - ieuj + 0TO,i(coswr in) i - isinwri„,i)) + l]x 
[—iu)(F 2 - ieuj + g in ,2{ COSUJT in,2 - «smwri„,2)) + 1] + 
w 2 g c ,ig c , 2 (cos(w(T c ,i + t C)2 )) - «sin(cj(r C! i + t c , 2 ))) = 

(11) 

which, after equating both the real and the imaginary 
parts with zero, provides for the implicit relations of uj 
and the system parameters 

-uj 2 P 1 P 2 + Q1Q2 = -^ 2 gc,\g c ,2 cos(cj(t Ci i + T C)2 )) 
0JP1Q2 +0JP2Q1 = uJ 2 g c ,igc : 2sin(uj(T Ct i +t c , 2 )), (12) 

where 

Pk = F k + g ln ^ k cos(uJT in . k ) 

Q k = eu 2 + g m , k uj sm(cJT in , fc ) - 1 (13) 

applies for k = 1,2. Squaring and adding the relations 
([HI) , one arrives at 



( W 2 P!P 2 - dQ,) 2 + w 2 (P!Q 2 + P2Q1) 2 = oj 4 g 2 c l9 2 



c,2i 

(14) 
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which can be used to express the cross-population cou- 
pling strengths in terms of w, while keeping the values for 
the subset of the intrinsic parameters <?in,fc) T in,fei &fc and 
Dk fixed. Obtained in a similar fashion, the analogous re- 
sult for the critical cross-population coupling delays may 
be written in the compact form 

Tc<1 + Tca = - arctanl^^^^^). (15) 

The last two equations combined define the curves in 
the appropriate delay-strength parameter plane. Bear in 
mind that Eq. (fT5"|) actually defines multiple branches of 
the Hopf bifurcation curves, these given by T Ci i+T Ci 2+i7r, 
where j — 0,1,2.... Naturally, the above relations further 
simplify once the inter-ensemble couplings are taken to be 
symmetrical and/or the populations' intrinsic parameters 
are assumed to be identical. Note that the expressions 
such as these could not be obtained if we were to retain 
the initial MF model ([7]) containing the full dynamics of 
the second order cumulants. As for the type of bifurca- 
tions whose location is indicated by (fT5| . the very form 
of the solution adopted for the characteristic equation 
is consistent with the Hopf bifurcations, though a rigor- 
ous proof would require one to verify whether the condi- 
tions on non-hyperbolicity, transversality and genericity 
are satisfied [12, 0, HH . Without entering into unneces- 
sary details, it suffices to say that the system ([5]) admits 
both the supercritical and subcritical Hopf bifurcations, 
whereby the former (latter) result in the creation of a 
stable (unstable) limit cycle. In addition, recall that ei- 
ther of these types can be direct or inverse [13] , depend- 
ing on whether an unstable two-dimensional manifold for 
the fixed point © appears or vanishes when crossing the 
bifurcation curve, respectively, having the fixed point un- 
fold on the unstable or the stable side. The results de- 
rived analytically are corroborated numerically by means 
of the DDE-biftool [49] , an adaptable package of Matlab 
routines suitable for handling the sets of DDE with con- 
stant delays. 



III. QUALITATIVE COMPARISON BETWEEN 
THE DYNAMICS OF THE EXACT AND THE 
APPROXIMATE SYSTEM 

For the systematic study, we first consider the layout 
with two populations made up of independent excitable 
elements {gi n ,i = ffm,2 = 0) subjected to a common, 
comparably small noise (Z?i = = 0.0001), whereby 
the cross-population coupling terms are taken to be sym- 
metrical, so one may introduce g Ct \ = g Ci i = g c and r Ci i = 
t c ,2 = t c . The parameters are such that for g c = 0, the 
populations exhibit the asymptotically (stochastically) 
stable equilibrium in the MF (exact) model. Though it 
appears marginal at first sight, the described setup is still 
important, since the MF model is here strongly indicated 
to match the behavior of the real system. In a sense, this 
scenario is reminiscent of a null-hypothesis, given that 



the stated parameters are fully compliant with the nom- 
inal conditions for the validity of the MF approximation. 
One would further expect to gain some insight into the 
phenomena occurring for the more complex system con- 
figurations, or may at least obtain a reference point to 
isolate the effects of certain parameters, such as gt n or 
Ti n . In the remainder, the bifurcation diagrams are ac- 
companied by the close-up views focused on the most 
relevant parameter domains, having those referred to in 
the text indicated by the representative symbols. Also, 
to distinguish between the different bifurcation curves, 
each is awarded with two types of indices. The +/ — 
sign specifies whether the curves coincide with the direct 
or inverse bifurcations, respectively, while the numerical 
index points to the order in which the given branches 
are encountered as the inter-population coupling delay is 
increased. 

Appreciating the bifurcation diagram in Fig. Ufa), a 
major point concerns the prediction on the existence of 
the critical strength go for the instantaneous couplings 
(t c = 0), where the stationary state loses stability. For 
g c < go, viz. the open circle in Fig. [TJb), the equilibrium 
is stable, whereas for g c > go (solid circle) there is only 
the oscillatory state. The bifurcation scenario coincides 
with the direct supercritical Hopf bifurcation, and the 
numerical simulations imply that the unstable manifold 
for the equilibrium m x ,i = m Xi 2 and m Vt i = m y ^ around 
.9c = .9o supports the oscillations in-phase, this being an 
example of synchronization between the units due to a 
common input. By the term "oscillations in-phase", it is 
meant that the MF approximation indicates a solution 
with the exact synchronization between the global vari- 
ables, which is stochastically perturbed in the exact sys- 
tem. What is described applies not only for r c = 0, but 
also holds in any instance when the curve ri,_ is crossed 
in the direction of increasing g c with r c kept fixed. How- 
ever, we note that there is an additional subtlety to this 
transition derived from an interplay with the fold-cycle 
bifurcation, a global event which cannot be accounted for 
by the present type of analysis. One finds that the sys- 
tem undergoes a global bifurcation around g c w 0.055, 
by which an unstable and a large stable limit cycle are 
born. This witnesses of an interval g c < go where the 
stationary and the oscillatory state coexist, with their 
attraction basins separated by the unstable limit cycle. 
Above go, the incipient limit cycle emerges around the 
former position of the equilibrium, but cannot fully grow 
with supercriticality, given that it is enclosed by the un- 
stable limit cycle created in the global bifurcation. At 
some point, the stable and the unstable limit cycle col- 
lide and disappear in an inverse fold-cycle bifurcation, 
which is indicative of a narrow g c interval around ti,_ 
corresponding to a bistable regime with a small and a 
large amplitude limit cycle, viz. the solid diamond in 
Fig. [Tfb). However, such bistability may be difficult to 
observe in the exact system for the sensitivity of the in- 
cipient cycle to stochastic perturbation, as even the very 
small noise amplitudes can prove sufficient to force the 
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ensuing orbits outside of its neighborhood. Away from 
<7o, the system's trajectory is eventually drawn to a dis- 
tant limit cycle attractor. 

To proceed with, we consider the effects of increasing 
t c under fixed coupling strength g c > g . Crossing the 
first bifurcation curve from below t c > ti—, viz. the 
domain indicated by an open triangle in Fig. QJb), the 
equilibrium is seen to regain stability via the inverse su- 
percritical Hopf bifurcation. Given the analogy of the 
underlying mechanisms, this could have been interpreted 
as a genuine example of the delay-induced amplitude 
death, if there were not for the large limit cycle which 
is unaffected by the local bifurcation. Instead, this sta- 
bility domain is characterized by the coexistence between 
the stationary and the oscillatory state. Nonetheless, en- 
hancing the delay above T2,+ gives rise to a region of 
instability, represented by the solid triangle in Fig. [TJb) , 
where one encounters only the two populations oscillat- 
ing in phase. Such an outcome is due to a supercritical 
Hopf bifurcation, which is reflected by the equilibrium 
gaining an unstable plane. Note that the analysis cannot 
extend to larger delays, since the underlying phenomena 
do not fall within the framework of the current study. 
It should be emphasized that the oscillation frequency 
of the MF model has been verified to match the one of 
the exact system almost perfectly. This point applies 
for two parameter domains highlighted by the solid cir- 
cle and the solid triangle in Fig. [IJb). Under Ti,-, the 
respective oscillation period of the approximate model is 
T ti MF = 3.664 in arbitrary units, whereas the associated 
average period for the exact system is T,_ex — 3.668. 
Likewise, in the domain instantiated by the solid trian- 
gle, T kMF = 3.874 and T k , EX = 3.869. The cited data 
witness that the MF model is able to predict the average 
frequency of macroscopic oscillations of the exact sys- 
tem with remarkable accuracy. From the aspect of com- 
parison between the real and the approximate systems, 
one should as well look back at the values of the critical 
strength go- The agreement here is weaker, whereby the 
MF model is found to overestimate the value. This is 
not unexpected, given that the local phenomena are me- 
diated by the background global bifurcation. Still, the 
tendency and rate by which go decreases with enhancing 
D is reflected reasonably well by the MF model. 

The main results in this Section concern the canon- 
ical setup involving two identical populations of inter- 
acting excitable neurons (g 

in.i — gin, 2 — 0-1), whereby 
the cross- pop ulation couplings are taken to be symmet- 
rical [3 |21| . The intrinsic ensemble parameters D = 
0.0001, Ti n = 0.3 warrant that the equilibrium is the only 
asymptotically (stochastically) stable state for the ap- 
proximate (exact) model. Inspecting the appropriate bi- 
furcation diagram in Fig. Ufa) , one readily realizes how, 
at variance with the previously discussed case, there is 
not one, but two scenarios for the destabilization of equi- 
librium. Which of the scenarios actually applies is con- 
tingent on the inter-population coupling strength g c : if 
9c < g'o, viz. Fig. [U[b), the equilibrium goes unstable 



via the direct supercritical Hopf bifurcation, while for 
g c > g' Q , the onset of the collective mode rests with the 
direct subcritical Hopf bifurcation. In the latter instance, 
where g c notably outweighs <?;„, an unstable limit cycle 
collapses at the fixed point making it unstable. Away 
from criticality, in the domain marked by the solid cir- 
cle in Fig(2fb) , the system's trajectory eventually gets 
drawn to a distant limit cycle attractor. Again, both 
the stable and the unstable limit cycle derive from the 
fold-cycle bifurcation, whereas the numerical simulations 
confirm that the unstable manifold of the equilibrium 
at (g c ,T c ) — (<?oj0) supports the symmetrical oscillatory 
state. Below the curve t^-, which is barely distinguish- 
able from the g c -axes in Fig. HtL), one finds a narrow 
interval of coupling strengths g c > go where the ema- 
nating branch of the unstable solutions apparently folds 
back. As a corollary, the system of coupled MF mod- 
els is seen to exhibit a bistable regime, such that the 
equilibrium and the collective mode coexist. However, 
such bistability is difficult to observe in the dynamics 
of the full system for the sensitivity of the equilibrium 
to stochastic perturbation. Interestingly, the approxima- 
tion for the critical coupling strength g' is significantly 
improved when compared to the previous system con- 
figuration, this possibly owing to the influence of the 
inter-ensemble interactions that were excluded earlier on. 
Crossing into the domain ti _ < t c < rg.-f represented 
by the open square in Fig. [2fb), the MF system under- 
goes an inverse subcritical Hopf bifurcation, such that 
the fixed point loses an unstable plane. Looking in a 
more general picture, this region of parameter space is 
supposed to be bistable between the equilibrium and 
the large limit cycle born via the global bifurcation. In 
parallel, the unstable limit cycle from the Hopf bifurca- 
tion should act like a threshold for switching between 
the two solutions. However, the stochastic component 
in the underlying dynamics prevents us from observing 
the bistable regime in the exact system. Above 72,+, the 
equilibrium loses stability, giving way to the limit cycle 
as the sole attractor of the system's dynamics. 

Next we turn to the sequence of bifurcations obtained 
for g c < g , which is a physically more plausible range 
since g c lies closer to gi n . Below T3 +, the stationary state 
is stable for both the real and the approximate system, 
with the appropriate parameter domain highlighted by 
the open up-triangle in Fig. HtL). Crossing T3, + from 
below, the system undergoes the supercritical Hopf bi- 
furcation, such that the equilibrium becomes unstable, 
and the emerging oscillations are symmetrical. An in- 
teresting point for the transition between the domains 
marked by the open and solid up-triangles in Fig. [2fb) is 
that for the moderate coupling strength, under not too 
large a noise the time lag turns out to be a necessary in- 
gredient should the equilibrium be destabilized. For the 
more comprehensive view, one again has to consider the 
effects of the interplay with the global fold-cycle bifurca- 
tion, whereby a general remark is that everything stated 
on the direct supercritical Hopf bifurcation regarding the 
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FIG. 1. (color online) The first few branches of the Hopf bifurcation curves r c (g c ) are shown for the MF based approximation, 
derived for the system of two symmetrically coupled populations made up of independent excitable units (gi„ — 0). The latter 
part is a deliberate idealization, intended to probe the predictions of the approximate system when all the conditions required 
for the validity of the MF model are exactly satisfied. The onset of the collective mode is an instance of synchronization among 
elements receiving the common input. The units within the ensembles are subjected to noise whose amplitude is D = 0.0001. 



diagram in Fig. [TJb) can carry over to this case. In brief, 
apart for the equilibrium, the system's phase space be- 
low t 3j+ also exhibits an unstable limit cycle enclosing 
the fixed point and a large stable limit cycle. Above the 
latter curve, the incipient limit cycle grows only until col- 
liding with the unstable one, both being annihilated in 
an inverse fold-cycle bifurcation. Then, all the trajec- 
tories are eventually drawn to the large limit cycle, left 
as the sole attractor. As for the predictions of the ap- 
proximate system, one stresses that there is an excellent 
agreement between the oscillating waveforms, in partic- 
ular when comparing the anticipated frequency with the 
average one for the real system, viz. T kt MF = 3.836 vs. 
T^.,ex = 3.833. This is illustrated in Fig. [3l showing 
side-by-side the sequences from the time series m Xj i(t) 
and Xi(t) for i = 1,2 below (top row) and above (bot- 
tom row) the curve T3, + . 

Further enhancing r c to step into the domain high- 
lighted by an open down-triangle in Fig. [2{b), one en- 
counters the bistable dynamics, such that the system, de- 
pending on the initial conditions, may display either the 
stationary or the oscillatory state. The area is bounded 



by T4 : _ from below and t$_ + from above. The found bista- 
bility regime is the consequence of the inverse subcritical 
Hopf bifurcation, where the emanating unstable cycle ef- 
fectively acts to stabilize the fixed point, allowing for it to 
coexist with the collective mode, the latter present due 
to the global bifurcation. The possibility of observing 
bistability in the exact system is likely facilitated by the 
unstable limit cycle, whose amplitude is sufficient to sep- 
arate more clearly between the attraction basins of the 
oscillatory solution and the equilibrium in spite of the 
stochastic perturbations induced by the relatively small, 
but non-negligible noise. The bistable regime is illus- 
trated in Fig. SI which demonstrates the coexistence of 
the stationary (top row) and oscillatory states (bottom 
row) for both the exact model and the MF approxima- 
tion. Note that the change in oscillating frequency in the 
real system, associated with crossing T4,_ from below, 
is well matched by the approximate system. Stepping 
into the domain rs. + < t c < Tg,-, marked by the solid 
down-triangle in Fig. HJb), the key change consists in 
the switch from the bistable to a monostable regime, the 
latter attributed the oscillatory state with the synchro- 
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FIG. 2. (color online) Bifurcation diagram in case of the symmetrical nonlinear cross-population coupling, as determined 
from the local bifurcation analysis of the approximate system made up of two identical interacting MF models. Unlike the 
conditions applied for Fig. [1] the excitable units within each ensemble of the corresponding exact system are connected in the 
all-to-all fashion, with the couplings characterized by the non-negligible strength and time lag. (a) and (b) show the T c (g c ) 
Hopf bifurcation curves, whereby the close-up view in (b) is focused on the most relevant region of the parameter plane. The 
intrinsic parameters D — 0.0001, Ti n = 0.3 and g on — 0.1 warrant that the isolated populations (single MF systems) exhibit the 
stochastically (asymptotically) stable stationary state. 



nization in-phase. The switch occurs as the system un- 
dergoes the direct supercritical Hopf bifurcation, which 
adds unstable directions, altering the stability of the fixed 
point. The change from the bistable to the monostablc 
regime occurs in the same fashion for the MF and the 
exact system. Setting r c above T6,-, see the domain rep- 
resented by the open diamond in Fig. [H[b), one finds 
the bistability regime reinstated. However, the transi- 
tion is accompanied by the modulation of the oscillating 
frequency, the point well reflected by the approximate 
system, viz. T 0j ex — 4.097 against T<>,MF = 4.119. In 
general, the increase of coupling delay is biased toward 
reducing the oscillating frequency. 

Note that the qualitatively similar sequence of bifur- 
cations is verified to persist in a range of values, if D 
and Ti n are set so to admit the stable stationary state as 
the sole attractor for the isolated populations. Nonethe- 
less, in order for this framework to reflect accurately the 
behavior of the exact system, one should not consider too 
large noise amplitudes. The perturbation from larger D 
may be envisioned as if leading to an effective broaden- 



ing of the bifurcation curves for the real system, which 
renders the entire sequence smeared, and the underlying 
qualitative changes difficult to discern. 
Asymmetrical coupling 

A question that naturally arises is whether and how is 
the physical picture so far modified by taking the asym- 
metrical, rather than the symmetrical cross-population 
coupling terms. We have examined two different sce- 
narios: by one, the couplings in cither direction retain a 
common time lag, but attain different strengths, whereas 
in the other, strengths are the same, but the transmis- 
sion delays are disparate. In the former case, the coupling 
strength in one direction, say g Ci \ is kept fixed, while g Ci i 
varies continuously. The bifurcation diagram in the t c - 
9c,2 plane is plotted in Fig. [SJa), whereby the intrinsic 
population parameters are identical to those stated in the 
caption of Fig. [5] One may immediately raise the issue 
of why is the bifurcation sequence profile much simpler 
compared to that in Fig. (5Ja) . The possible reason lies in 
that for the cross-population couplings asymmetrical by 
strength, the system's behavior is predominantly influ- 
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FIG. 3. The exact (left column) and the approximate sys- 
tem (right column) are demonstrated to undergo the direct 
supercritical Hopf bifurcation when crossing the curve ti,+ 
from Fig. HJb). (a) and (b) show that below the curve 
(g c = 0.16, t c = 0.06), the fixed point is stochastically stable 
for the exact, and asymptotically stable for the approximate 
system, respectively. The onset of oscillations above the curve 
(g c = 0.16, t c = 0.14) is illustrated for the exact system in (c), 
and the approximate system in (d). The intrinsic population 
parameters are set to D = 0.0001, gi n = 0.1 and T in = 0.3. 
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FIG. 4. Illustration of how the same bistable regime, char- 
acterized by coexistence between the stationary and the os- 
cillatory state, is exhibited both by the exact (left column) 
and the approximate system (right column). The top row 
indicates the corresponding stochastically and asymptotically 
stable fixed point, whereas the bottom row shows the two pop- 
ulations oscillating in-phase. The coupling strength and delay 
(g c , t c ) = (0.14, 0.22) lie within the domain highlighted by the 
open down-triangle in Fig. [2jb). The values for the intrinsic 
parameter subset are D — 0.0001, gm = 0.1 and T in = 0.3. 



enced by the global bifurcation phenomena dependent on 
g Ci i and g c ,2- Nonetheless, one cannot neglect some qual- 
itative resemblance between the dynamics of the MF and 
the exact system. For instance, below in Fig. Ufa), 
the equilibrium is stable for either system, but partici- 
pates in the bistable regime. Along with the stationary 
state, one also finds an oscillatory state where the two 
populations are locked with a constant phase shift. This 
collective mode can only be attributed to the global bi- 
furcation events. Crossing ri. + from below results in the 
creation of a limit cycle, leaving the equilibrium unsta- 
ble. Both the real and the approximate model exhibit a 
single attractor supporting the phase-locked oscillations 
between the two populations, whereby the underlying fre- 
quencies are well matched, viz. T,^mf — 4.281 against 
T, : ex = 4.302. Notably, the oscillation waveforms above 
Ti + are more complex than those below, and bear the 
initial signatures of the quasiperiodic behavior. It has 
to be stressed that the qualitative resemblance between 
the dynamics of the exact and the approximate system 
heavily depends on how close is g Ct \ to <?j„. In Fig. 03a), 
g c i — 0.05 is comparably small to g^ n = 0.1. Should g c> \ 
approach g^ n or exceed it, the effects of the global bi- 
furcation phenomena become overwhelming, spoiling the 
predictions made by MF-based approximation. 

We also briefly touch upon the setup where the cross- 
population couplings exhibit the disparate time lags, but 
attain the same coupling strength. Again, all the internal 
population parameters are equal to those linked to Fig. 



[2j whereas the notation on the asymmetrical coupling 
parameters is analogous to that used in the previous lay- 
out. The appropriate bifurcation diagram in the T c .2-g c 
plane is displayed in Fig. EJb). Compared to Fig. EJa), 
we learn how the main difference between this case of 
asymmetrical couplings and the case with symmetrical 
interaction lies in the domain of small delays. In partic- 
ular, the destabilization of equilibrium occurs solely via 
the supercritical Hopf bifurcation, whereas the scenario 
involving the subcritical Hopf bifurcation is absent. This 
picture seems to be independent on the relation between 
the fixed time lag t c ,i and Ti n . 

Though it is not within the scope of the current study, 
one should still mention that the methods discussed 
can also be implemented for the scenarios where the 
two populations exhibit different types of kinetics, e.g. 
if one is made up of excitable, and the other of self- 
oscillating units. In this scenario, one effectively exam- 
ines the interaction between the noise-induced and the 
noise-perturbed oscillations. The corresponding bifurca- 
tion diagram is not too distinct from the one in Fig. [S[b), 
except that the pattern of bifurcation curves is less dense. 
The critical coupling strength analogous to go is naturally 
smaller than the one for the interacting excitable popu- 
lations. Nevertheless, this setup is distinguished from 
those considered earlier in that the unstable manifold 
of the equilibrium supports the onset of the collective 
mode with the phase-locked rather than the in-phase os- 
cillations, such that the firing of the ensemble with self- 
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FIG. 5. (color online) Results of the local bifurcation analysis carried out on the approximate system for the two cases of 
asymmetrical cross-population couplings, are presented in the delay-strength parameter plane, (a) refers to the setup with the 
disparate coupling strengths, holding g c ,i — 0.05 and letting g Ct 2 vary continuously, (b) is obtained for the uneven time lags, 
with Tc,i = 0.6 fixed and t Cj 2 allowed to change. The intrinsic parameters D = 0.0001, Ti n — 0.3, gin = 0.1 are identical for 
both populations and warrant that the corresponding isolated system would exhibit the stationary state. 



oscillating neurons precipitates the firing of the ensemble 
containing the excitable neurons. 

IV. SUMMARY AND DISCUSSION 

In the present paper, we have pursued the analysis of 
the MF based approximation intended to accurately re- 
flect the macroscopic behavior of two delay-coupled pop- 
ulations of stochastic excitable units in terms of the sta- 
bility of the stationary state, the scenarios for the onset 
and the suppression of the collective mode, as well as 
the possibility of admitting bistable regimes, where the 
equilibrium and the oscillatory state are found to coex- 
ist. The described layout deserves attention, since it can 
be interpreted as the minimal model for the "network of 
networks" , the configuration often brought into context 
of biological systems whose function relies on generation 
and adjustment between the multiple collective rhythms. 
The important ingredients of the exact system we con- 
sider include two types of delayed interactions, whereby 
those within the ensembles are assumed to be linear, and 



the inter-ensemble ones, mediated by the appropriate 
global variables, are taken to be nonlinear. The corre- 
sponding approximate system is built by coupling the 
two MF models, derived to describe the activity of single 
populations. Such a framework follows the general idea 
that any ensemble of oscillating units exhibiting the col- 
lective mode can be treated as the macroscopic oscillator. 
The MF model integrates the cumulant approach with 
the Gaussian approximation, whereby the latter holds 
exactly if three conditions are satisfied. These include 
the thermodynamic limit N — > oo regarding the ensem- 
ble size, the negligible noise amplitude D << 1, as well 
as the negligible interaction between the units < < 1 . 
However, it turns out that the approximate system is still 
able to predict with sufficient accuracy the behavior of 
large, but finite populations with non-negligible internal 
interactions, provided that the natural requirement for 
not too large a noise amplitude is met. 

By stating the results in broad terms, the intention 
has been to stress their applicability to the class of sys- 
tems made up of type II excitable units. Nonetheless, one 
recognizes that valuable motivation for the study comes 
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from the field of neuroscience, which goes beyond the 
adopted model of local dynamics or the fashion in which 
the interactions are introduced. The methods for pro- 
viding the reduced descriptions of the behavior of large 
neural assemblies are typically cast in the categories of 
the neural-mass and the MF models, whereby the former 
neglect, and the latter take into account the distribution 
of individual neuron states over the ensemble. In these 
terms, the model considered here interpolates between 
the two classes. Recall that we have introduced an ad- 
ditional approximation on the second order moments to 
translate the original MF system ([7]) into the form incor- 
porated in ([8]), with the latter preferred as it allows for 
the analytical tractability of the subsequent local bifur- 
cation analysis. 

An inference from such an analysis is that the approxi- 
mate system can undergo direct and inverse supercritical 
or subcritical Hopf bifurcations, such that the direct (in- 
verse) ones lead to the destabilization (stabilization) of 
the stationary state. The complex bifurcation sequence 
under variation of cross-population coupling strengths 
and delays is found to depend on the details of the sys- 
tem configuration, like the symmetrical or asymmetrical 
character of the bidirectional interaction between the en- 
sembles. The main set of results refers to the symmetrical 
case, where it is demonstrated that the equilibrium may 
lose stability according to two different scenarios. One in- 
volves a direct supercritical Hopf bifurcation and can be 
achieved for instantaneous couplings solely by increasing 
<7 C , whereas the other scenario unfolds via the direct sub- 
critical Hopf bifurcation. The latter involves an interest- 
ing point that for strengths g c ~ gi n one finds a time-lag 
threshold necessary to destabilize the equilibrium. In- 
creasing t c , there are parameter domains bounded from 
below (above) by the curves indicating subcritical (su- 
percritical) bifurcations, where the stability of station- 
ary state is regained. In many of such instances, the 
system is actually bistable, exhibiting coexistence be- 
tween the equilibrium and the oscillatory state. This is 
a corollary of an interplay with the global fold-cycle bi- 



furcation, as the large stable limit cycle born in this way 
remains unaffected by the local phenomena. Note that 
the global events may influence the system dynamics in 
several other instances. In particular, an unstable limit 
cycle created in a fold-cycle bifurcation may destabilize 
the fixed point in a direct subcritical Hopf bifurcation or 
may limit the growth of an incipient limit cycle following 
the direct supercritical Hopf bifurcation. By numerical 
simulation, we have verified that the parameter domains 
of stability or instability of equilibrium for the exact sys- 
tem are reproduced by the approximate one with high 
accuracy. In addition, it has been shown that the av- 
erage oscillation frequency for the global variable of the 
exact system is well matched by that of the correspond- 
ing MF variable. In the exact system, the ability to ob- 
serve the bistable regimes, where the unstable limit cycle 
act as a threshold between the equilibrium and the large 
cycle, is contingent on the noise amplitude. In general, 
the predictions of the approximate system are better if 
the deterministic component, governed by the coupling 
strengths and delays, prevails over the stochastic compo- 
nent in the dynamics of the exact system. An interesting 
study complementary to the present one would be to ex- 
amine whether the MF based model may reproduce the 
forms of synchronization between the generated collec- 
tive rhythms the way they are exhibited by the exact sys- 
tem. These could include the in-phase and antiphase syn- 
chronization or the phase-locked states, as well as their 
coexistence. The preliminary results implementing the 
H-function approach suggest that the approximate sys- 
tem may account for the stability of the synchronization 
regimes and provide indications on the possible multista- 
bility. 

ACKNOWLEDGMENTS 

This work was supported in part by the Ministry of 
Education and Science of the Republic of Serbia, under 
project Nos. 171017 and 171015. 



[1] D. Golomb, D. Hansel, and G. Mato, in Neuro- 
informatics and Neural Modeling, edited by F. Moss and 
S. Gielen, Handbook of Biological Physics, Vol. 4 (El- 
sevier, Amsterdam, 2001); J. L. P. Velazquez and R. 
Wennberg eds., Coordinated Activity in the Brain: Mea- 
surements and Relevance to Brain Function and Behav- 
ior, (Springer, New York, 2009); G. Buzsaki, Rhythms of 
the Brain (Oxford University Press, Oxford, 2006). 

[2] P. A. Tass, Phase Resetting in Medicine and Biol- 
ogy: Stochastic Modelling and Data Analysis, (Springer, 
Berlin Heidelberg, 2007). 

[3] S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. 
Yagita, M. Kobayashi, and H. Okamura, Science 302, 
1408 (2003). 

[4] L. Glass and M. C. Mackey, From Clocks to Chaos: The 
Rhythms of Life (Princeton University Press, Princeton, 



NJ, 1988). 

[5] S. H. Strogatz, D. M. Abrams, A. McRobie, B. Eckhardt, 
and E. Ott, Nature (London) 438, 43 (2005); Z. Neda, 
E. Ravasz, Y. Brechet, T. Vicsek, and A.-L. Barabasi, 
Nature (London) 403, 849 (2000). 

[6] I. Kiss, Y. Zhai, and J. Hudson, Science 296, 1676 (2002). 

[7] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Phys. Rev. 
Lett. 76, 404 (1996). 

[8] M. Y. Kim, R. Roy, J. L. Aron, T. W. Carr, and I. B. 
Schwartz, Phys. Rev. Lett. 94, 088101 (2005). 

[9] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares, and 
C. S. Zhou, Phys. Rep. 366, 1 (2002); A. Pikovsky, M. 
G. Rosenblum, and J. Kurths, Synchronization, A Uni- 
versal Concept in Nonlinear Sciences (Cambridge Uni- 
versity Press, Cambridge, 2001). 
[10] N. Buric and D. Todorovic, Phys. Rev. E 67, 066222 



14 



(2003) ; M. Dhamala, V. K. Jirsa, and M. Ding, Phys. 
Rev. Lett 92, 074104 (2004); O. V. Popovych, C. Haupt- 
mann, and P. A. Tass, Phys. Rev. Lett. 94, 164102 
(2005); M. D. McDonnell and L. M. Ward, Nat. Rev. 
Neurosci. 12, 415 (2011). 

[11] D. V. Ramana Reddy, A. Sen, and G. L. Johnson, Phys. 
Rev. Lett. 80, 5109 (1998); S. H. Strogatz, Nature 394, 
316 (1998).. 

[12] A. Pikovsky, A. Zaikin, and M. A. de la Casa, Phys. Rev. 
Lett. 88, 050601 (2002); B. Lindner, J. Garcia-Ojalvo, A. 
Neiman, and L. Schimansky-Geier, Phys. Rep. 392, 321 

(2004) ; N. Buric, K. Todorovic, and N. Vasovic, Phys. 
Rev. E 75, 026209 (2007); Q. Wang, G. Chen, and M. 
Perc, PLoS ONE 6, el5851 (2011). 

[13] Y. Baibolatov, M. Rosenblum, Z. Z. Zhanabaev, M. Kyz- 
garina, and A. Pikovsky, Phys. Rev. E 80, 046211 (2009). 

[14] Y. Kawamura, H. Nakao, and Y. Kuramoto, Phys. Rev. 
E 75, 036209 (2007); Y. Kawamura, H. Nakao, K. Arai, 
H. Kori, and Y. Kuramoto, Phys. Rev. Lett. 101, 024101 
(2008); H. Kori, Y. Kawamura, H. Nakao, K. Arai, and 
Y. Kuramoto, Phys. Rev. E 80, 036207 (2009). 

[15] Y. Kawamura, H. Nakao, K. Arai, H. Kori, and Y. Ku- 
ramoto, CHAOS 20, 043109 (2010). 

[16] E. Ott and T. M. Antonsen, CHAOS 18, 037113 (2008). 

[17] E. Barreto, B. Hunt, E. Ott, and P. So, Phys. Rev. E 77, 
036107 (2008). 

[18] E. Montbri, J. Kurths, and B. Blasius, Phys. Rev. E 70, 
056125 (2004). 

[19] M. G. Rosenblum and A. S. Pikovsky, Phys. Rev. Lett. 

92, 114102 (2004). 
[20] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. 

Wiley, Phys. Rev. Lett. 101, 084103 (2008). 
[21] S. Olmi, A. Politi and A. Torcini, EPL 92, 60007 (2010). 
[22] E. M. Izhikevich, Dynamical Systems in Neuroscience: 

The Geometry of Excitability and Bursting (MIT Press, 

Cambridge, Massachusetts, 2007). 
[23] P. S. Skardal and J. G. Restrepo, Phys. Rev. E 85, 016208 

(2012). 

[24] H. Hasegawa, Physica D 237, 137 (2008). 

[25] B. Lindner, J. Garcia-Ojalvo, A. Neiman, and L. 

Schimansky-Geier, Phys. Rep. 392, 321 (2004). 
[26] R. A. Stefanescu and V. K. Jirsa, PLoS Comput Biol. 4, 

el000219 (2008). 
[27] O. David and K. Friston, Neuroimage 20, 1743 (2003). 
[28] R. Rodriguez and H.C. Tuckwell, Phys.Rev.E, 54, 5585 

(1996). 

[29] S. Tanabe and K. Pakdaman, Phys.Rev. E, 63, 031911 
(2001) . 

[30] M.A. Zaks, X. Sailer, L. Schimansky-Galer and A.B. 
Neiman, Chaos, 15,026117 (2005) . 



[31] M. A. Zaks, A. B. Neiman, S. Feistel, and L. Schimansky- 
Geier, Phys. Rev. E 68, 066206 (2003). 

[32] H. Hasegawa, Phys.Rev. E, 68, 041909 (2003). 

[33] H.Hasegawa, Phys.Rev. E, 72 021912 (2004). 

[34] H.Hasegawa, Phys.Rev. E, 72, 021911 (2004). 

[35] N. Buric, D. Rankovic, K. Todorovic, and N. Vasovic, 
Physica A 389, 3956 (2010). 

[36] N. Buric, K. Todorovic, and N. Vasovic, Phys. Rev. E 
82, 037201 (2010). 

[37] F. M. Atay, editor, Complex Time-Delay Systems: 
Theory and Applications, (Springer, Berlin Heidelberg, 
2010); E. Scholl, G. Hiller, P. Hovel and M. A. Dahlem, 
Phil. Trans. R. Soc. A 367, 1079 (2009). 

[38] S. Luccioli and A. Politi, Phys. Rev. Lett. 105, 158104 
(2010); M. Rosenblum, A. Pikovsky, Phys. Rev. E 70, 
041904 (2004). 

[39] R. E. Lee DeVille, E. Vanden-Eijnden and C. B. Muratov, 
Phys. Rev. E 72, 031105 (2005); C. B. Muratov and E. 
Vanden-Eijnden, CHAOS 18, 015111 (2008). 

[40] S. Tanabe, K. Pakdaman, Phys. Rev. E 63, 031911 
(2001). 

[41] M. A. Zaks, X. Sailer, L. Schimansky-Galer, and A. B. 
Neiman, Chaos 15, 026117 (2005). 

[42] L. Arnold, Random Dynamical Systems (Springer- 
Verlag, Berlin, 1998). 

[43] C. W. Gardiner, Handbook of Stochastic Methods for 
Physics, Chemistry and the Natural Sciences, (Springer- 
Verlag, Berlin, 1985). 

[44] I. Franovic, K. Todorovic, N. Vasovic and N. Buric, Phys. 
Rev. Lett. 108, 094101 (2012). 

[45] J. Hale and S. V. Lunel, Introduction to Functional Dif- 
ferential Equations ( Springer- Verlag, New York, 1993). 

[46] S. A. Campbell in Handbook of Brain Connectivity, edited 
by R. Mcintosh and V. K. Jirsa. (Springer- Verlag, Berlin 
Heidelberg, 2007); S. A. Campbell in Delay Differential 
Equations: Recent Advances and New Directions, edited 
by B. Balachandran, T. Kalmr-Nagy and D. Gilsinn 
(Springer- Verlag, New York, 2009). 

[47] S. Wiggins, Introduction to Applied Nonlinear Dynam- 
ical Systems and Chaos, 2nd ed. (Springer, New York, 
Cambridge, 2000). 

[48] Y. A. Kuznctsov, Elements of the Applied Bifurcation 
Theory, 3rd ed. (Springer- Verlag, New York, 2004). 

[49] K. Engelborghs, T. Luzyanina, and G. Samaey, Technical 
Report TW-330, Department of Computer Science, K. 
U. Leuven, Leuven, Belgium (2001); K. Engelborghs, T. 
Luzyanina and D. Roose, ACM Trans. Math. Softw. 28, 
1 (2002). 



